Manifold Learning Approach for Chaos in the Dripping Faucet 
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Dripping water from a faucet is a typical example exhibiting rich nonlinear phenomena. For such a system, the 
time stamps at which water drops separate from the faucet can be directly observed in real experiments, and the 
time series of intervals r„ between drop separations becomes a subject of analysis. Even if the mass m„ of a drop 
at the onset of the n-th separation, which cannot be observed directly, exhibits perfectly deterministic dynamics, 
it sometimes fails to obtain important information from time series of t„. This is because the return plot t„_i vs. 
t„ may become a multi-valued function, i.e., not a deterministic dynamical system. In this paper, we propose a 
method to construct a nonlinear coordinate which provides a "surrogate" of the internal state m n from the time 
series of t„. Here, a key of the proposed approach is to use ISOMAP, which is a well-known method of manifold 
learning. We first apply it to the time series of t„ generated from the numerical simulation of a phenomenological 
mass-spring model for the dripping faucet system. It is shown that a clear one-dimensional map is obtained by 
the proposed approach, whose characteristic quantities such as the Lyapunov exponent, the topological entropy, 
and the time correlation function coincide with the original dripping faucet system. Furthermore, we also 
analyze data obtained from real dripping faucet experiments which also provides promising results. 

PACS numbers: 05.45.-a, 05.45.Tp, 05.10.-a 



I. INTRODUCTION 

Dripping of water from a faucet is ordinarily seen in our 
daily life. At first glance, such a motion of dripping looks 
very common. It provides, however, a variety of rich non- 
linear dynamics including the period-doubling bifurcation to 
chaos, intermittency, crisis, hysteresis, and etc. In particular, 
Robert Shaw and his collaborators [1] first found that there is 
a clear transition from a periodic motion to low-dimensional 
chaos by investigating the time intervals r„ between dripping 
separations from the faucet, both theoretically and experimen- 
tally. 

The dripping water is fluid dynamics, i.e., ideally described 
as an infinite-dimensional dynamical system. But as far as 
the dynamics is confined within a low dimensional attrac- 
tor, it can be modeled by a class of phenomenological mod- 
els called "mass-spring" systems Since the pioneering 
work of Shaw et al., many versions of the mass-spring system 
for the dripping faucet have been proposed. Among them, 
Kiyono and Fuchikami [2] significantly improved the mass- 
spring system on the basis of both, numerical simulations of 
fluid dynamics JH and real experiments fl. They showed that 
their model can systematically explain various aspects of the 
complex behaviors observed in the real dripping faucet exper- 
iments. 

One of the most prominent aspects of the Kiyono- 
Fuchikami model is that the essential feature of chaos in the 
dripping faucet is exactly represented as a one-dimensional 
map. More precisely, the mass m„ at the moment of the n-th 
separation of a drop from the faucet obeys a one-dimensional 
mapping dynamical system, i.e., there exists a deterministic 
scalar function /(•) such that m n = /(m„_i). 



In general, however, not all state variables are observable 
in real experiments. In the case of the dripping faucet system, 
for example, it is very difficult to observe the mass m n of a 
drop in a direct way. Instead, time intervals r„ between suc- 
cessive drop separations can be recorded in real experiments. 
As investigated by Shaw et al., depending on the degree of 
flux of water, the return plot r„_i vs. r„ also shows a clear 
functional relationship. At the same time, however, they have 
also shown that it often takes the form of a multi-valued rela- 
tion. Namely, there are two or more candidates of t„ against 
a single value of r„_i, which prevents us from interpreting 
the dripping faucet as a simple one-dimensional mapping sys- 
tem. This multi-valuedness problem often occurs in general 
chaotic dynamical systems such as the Kuramoto-Sivashinsky 
equation 150 . 

On the other hand, the existence of the one-dimensional 
map / associated with the mass m„ means that an embedding 
of the dripping-time interval r„ into a sufficiently high, say d- 
dimensional Euclidean space R d as s n = (T„_<f+i, ...,t„_i,t„), 
is lying on a one-dimensional manifold S dE. d . Then, a point 
s„ E S obeys a deterministic law as s„ = F(s„_i) where F(-) 
is a (i-dimensional vector valued function whereas the rela- 
tionship between r„_i and t„ is a multi-valued one. Actu- 
ally, in the case of the Kiyono and Fuchikami's mass-spring 
model, embedding t„ into a three-dimensional space gener- 
ally results a filament-like one-dimensional manifold without 
crossing. Therefore, if a new coordinate u is spanned along 
T, which plays the role of a surrogate variable for the inter- 
nal state m n , then we obtain a more simplified expression as 
m„ = g(M„_i) where g(-) is a scalar function of u. Even if the 
original one-dimensional mapping system m„ = /(w„_i) is 
not available, important dynamical features can be obtained 



from the mapping associated with the surrogate variable u as 

u„ = g(u„_i). 

To identify lower dimensional representations of the dy- 
namics, dimension reduction methods can be employed. Di- 
mension reduction is an important task of data (pre-) process- 
ing with applications in pattern and speech recognition, image 
processing, bioinformatics, psychology, etc. Linear subspaces 
containing or approximating the available data can be identi- 
fied using Principal Component Analysis (PCA), Independent 
Component Analysis (ICA), and many useful methods in vari- 
ous fields |6[] . However, when the data set of interest is located 
on or close to a (sub-) manifold with significant curvature, the 
applicability of these linear methods is limited and nonlinear 
dimension reduction methods have to be employed. 

Recently, in the field of statistical machine learning, meth- 
ods of manifold learning have been developed for providing a 
low-dimensional representation when data is lying on a non- 
linear low dimensional manifold embedded in a high dimen- 
sional Euclidean space. A number of methods have been pro- 
posed, and in the present study, we employ ISOMAP for 
such a purpose. ISOMAP, which is an abbreviation of the term 
"isometric feature mapping", is a method of manifold learning 
where the geodesies between training samples are employed 
as the dissimilarity information in multi-dimensional scaling 
(MDS) H. 

In this paper, we demonstrate that ISOMAP is very useful 
to extract a surrogate state variable u and to construct a well 
defined one-dimensional map g( -) for both, numerical and real 
experimental data. It is shown that dynamical characteristics 
such as the Lyapunov exponent and the time correlation func- 
tion can be computed from g(-). 

The present paper is organized as follows. In Section II, 
we explain the dripping faucet system. In Section III, we first 
introduce the method of ISOMAP, then we apply it to data 
generated from the mass-spring model mentioned in the pre- 
vious section. Finally, in Section IV, we give a summary and 
discuss possible directions of future research. 
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FIG. 1 . (a) Snapshots of a high speed movie showing the separation 
of a drop, (b) Illustration of the mass spring model. 



oscillations. 



B. Mass-spring model 



II. DRIPPING FAUCET SYSTEM: MODEL AND 
EXPERIMENT 

A. Basic mechanism 

Let us begin with a brief introduction of the basic mecha- 
nism how a water drop separates from a faucet. Figure Q~|(a) 
shows a snapshot at just the moment when a water drop sep- 
arates from a burette in experiments. Here, the shape of the 
drop is determined by the balance between the surface tension 
and the weight of water. When increasing the mass of a drop 
by injecting water, the following processes are repeated with 
time, (i) A "neck" which connects between the drop and the 
faucet is formed by the break of the balance between the ten- 
sion and the mass of water, (ii) When the weight reaches a 
critical value, the neck is broken, i.e., a portion of the drop 
separates from the faucet, (iii) Just after its separation, the re- 
mainder of the drop rapidly shrinks by the surface tension to 
the upward direction, (iv) Finally, the drop grows again with 



Based on observations as mentioned in the previous sub- 
section, the following equations of motion can be considered 
as a phenomenological model for the dripping faucet experi- 
ment in 

m'x + xm — —kx — yx + mg, (1) 
m = Q (const.). (2) 

Such a model is called a mass-spring model and its schematic 
illustration is depicted in Fig. [T](b). Here, x is the vertical po- 
sition of the forming drop to the downward direction, m is its 
mass, g is the gravitational acceleration, k is the stiffness of the 
spring, and y is the damping parameter. The restoring force 
given by the surface tension is represented as a spring force in 
Eq. (03, and the mass of the forming drop linearly depends on 
time with the rate Q as described in Eq. ©, because a drop 
grows with time due to the influx of water from the faucet. It 
is also assumed that when the position of the drop reaches the 
critical point x c , the drop loses its mass by Am due to the sep- 
aration of a portion of the drop and this portion falls into the 
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ground. In spite of its simplicity, this model can explain many 
dynamical aspects of the real dripping faucet QJJ] . 

Kiyono and Fuchikami improved the above phenomenolog- 
ical model [2] based on the knowledge of the numerical simu- 
lations of fluid dynamicsi [3] and real experiments of the drip- 
ping faucet [4]. They first modified the equation of motion in 
Eq. (Q} as 



mx + (i - Vo)ih — —kx — yx + mg, 



(3) 



where vo is the velocity of the influx of water. Note here that 
there is a relation between Vo and Q in Eq. as Q = na 2 vo 
where a is the radius of the faucet. Then, based on their real 
experiments, they considered that the stiffness k in Eq. (f3]) also 
depends on the mass of the drop as: 



k(m) 



-11.4m + 52.5 

(m > m c ), 



(m < m c ), 
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where m c = 4.61. Equation (0]l means that when the 
mass m amounts to m c , the value of the stiffness be- 
comes zero, then the drop undergoes free-fall. In their 
experiments, the units of the length, time, and mass are 
chosen as (y/pg) 1/2 (=0.27cm), (F/pg 3 ) 1/4 (=0.017 sec) and 
p(y/pg) 3 l 2 ={0.020g), respectively, where F is the surface ten- 
sion and p is the density. Using these units, parameters are 
set to y = 0.05, g = 1, x c = 5.5, Am = 0.8w - 0.3, and 
a = 0.916, and the constants in Eq. (01 are also determined 
from their experiments. They also assumed that just after a 
portion of the drop separates, the position and velocity are re- 
set to x — xq — 2.0, x — io = 0. 

Figure [2] shows a trajectory of the above mentioned model 
with vo = 0.1 130 after some transient. In Fig. |2](a) we can see 
that the trajectory is tracing a chaotic attractor. In real experi- 
ments, however, it is in general impossible to observe all state 
variables of the system. In the case of the dripping faucet, time 
intervals r„, (n — 1, 2, ...) between successive drop separations 
are observed in experiments. How to determine r„ from the 
signal of the position x(t) is depicted in Fig. |2](b). Here the 
variable m„ is the value of the mass at the moment of the n-th 
drop separation. 

Figure [3] (a) shows the return plot m„_i vs. m„ for vo 

0. 1 HO We can see that there is a clear scalar function be- 
tween m„_i and m„ as m„ = f(m n -i). In Figure [3] (b), how- 
ever, the return plot r„_i vs. r„ is a multi-valued function, 

1. e., the right-hand side of the return plot shows the 1 to 2 val- 
ues. From a different viewpoint, if we regard this return plot 
as the time-delay embedding of r„ into the two-dimensional 
plane R 2 as s n = (t„_i,t„) denoting the manifold on which 
the states s n are lying as S, then, we can see that there is a de- 
terministic relationship between s„_i and s„ as s„ = F(s„_i) 
in R 2 . 



C. Bifurcation Structure 

We also investigated how the statistical property of the 
mass-spring model (Eqs. <[3j and ©) depends on the water 



FIG. 2. (a) Chaotic trajectory of the Kiyono-Fuchikami model, (b) 
Time series of drop position x and the time interval t between drop 
separations. 
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FIG. 3. (a) Return plot of m„_i vs. m„. (b) Return plot of t„_i 
vs. t„. Movements of points as the time delay embedding vector 
s„ = (t„_i,t„) are also depicted. 



influx vo, and the result is shown in Fig.|4](a). One can see rep- 
etitions of the period doubling bifurcation route to chaos, as 
well as periodic windows and their reverses as increasing vq. 
We also made experiments to check whether the true dripping 
faucet system also exhibits this bifurcation structure. Figure 
|4](b) shows a time series of the time intervals of drop separa- 
tions over a long time period. Here, in our experiments, the 
surface of water of the bath decreases very slowly because no 
water is supplied from outside, which plays a role of chang- 
ing the water influx. Therefore, this figure represents a kind 
of "bifurcation" diagram. One can see that there is a signifi- 
cant qualitative similarity between the numerical simulations 
(Fig-H(a)) and the real experiments (Fig.|4](b)). 
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FIG. 4. (a) Bifurcation diagram of the Kiyono-Fuchikami model. 
Note that the unit of time is rescaled, e.g., r = 14 corresponds to 
0.238 sec. (b) Time series of a real dripping faucet experiment. 



III. EXTRACTING ONE-DIMENSIONAL MAPS OF 
INTERNAL STATE VARIABLES BY ISOMAP 



As shown in the previous section, this dripping model is es- 
sentially described by a one-dimensional map f(m). Its exper- 
imental observables r„ also show a one-dimensional filament 
in the cZ-dimensional space (t„_^ + i, T n -d+i, ■ ■ ■ , t„), but the re- 
lationship between r„_i vs. r„ is not always given by a one- 
dimensional map directly. In this paper, we discuss the case 
d — 2 as shown in Fig|3] This suggests that the dripping-time 
interval t isn't appropriate for the simple description of the 
dynamics. If we can construct a new coordinate u along the fil- 
ament, the dynamics must be described by a one-dimensional 
map u n = g(u n -\) and easily analyzed using the theory for 
one-dimensional maps. In this section, we try to construct a 
new coordinate u by applying ISOMAP to the time series r„ 
and get the one-dimensional map g(u). In addition, we test 
whether we can recover the statistical properties of the drip- 
ping faucet system from the one-dimensional map g(u). 



A. ISOMAP 

ISOMAP is one of several widely used low-dimensional 
embedding methods, which is an extension of classical MDS 
(multi-dimensional scaling) |7y. MDS seeks a low dimen- 
sional representation of the sample points. This is achieved 
by plotting data points in a low dimensional space preserv- 
ing the "dissimilarity" (generalized distance) between sample 
points (in the original higher dimensional space) as much as 
possible. For the dripping time series t„, when we use the Eu- 
clidian distances (pairwise distances between sample points) 
in the (t„_i,t„) plane as the dissimilarity, we shall not ob- 
tain a one-dimensional embedding, because the sample points 
are located on a curved filament. In ISOMAP, the geodesic 
distance on the low dimensional structure instead of the Eu- 
clidean distance is used and then we can embed the sample 
points into the one-dimensional space and get the new coordi- 
nate u. A concrete procedure is described as follows. 

First, we compute the geodesic distance djj between z'-th 
and j'-th sample points which is approximated with the short- 
est path from one to the other on the neighboring graph Q. 
This graph Q is constructed by locally connecting among sam- 
ple points (we employ the Euclidean distance to construct k- 
nearest neighbors). This procedure with k = 3 is illustrated 
in Fig. 13 a). The point i is connected to i\, iz, h, which are k- 
nearest neighbors of the point i and all sample points are con- 
nected in the same manner. The distance djj between points i 
and j is not defined by the Euclidean distance (dashed line), 
but by the shortest path distance (heavy solid line). As a result, 
dij is an approximation of the geodesic distance on the mani- 
fold. Figure|5jb) shows the graph Q of the dripping faucet (0 
and © for v = 0.1 130, AT = 500 and k = 4. And d u is ob- 
viously a good approximation of the geodesic distance along 
the filament. 

The bifurcation diagram for this dripping faucet is shown in 
Fig. [4] When the attractor has n-bands or exhibits strong inter- 
mi ttency, the neighboring graph @ may be separated into a few 
clusters and then we cannot estimate the geodesic distances 
dij. In this case we adopt the shortest connection among the 
clusters to construct the global graph Q. A example for two- 
band chaos (vo = 0.1128) is shown in Fig|5jc). 

The second step is MDS whose purpose is to place a set of 
new points m,, (z = 1, 2, . . . , N) in a low dimensional space 
so that the dissimilarities dij in the original state are well- 
approximated by \ui - Uj\, i.e. we find points m, that min- 
imize ^ (dij - |«j - Uj\j . In MDS, for the sample points 

y 

Si, (i= 1, 2, . . . , N) in the ra-dimensional Euclidean space, we 
define the square distance matrix D as Dij = |s, - Sj\ 2 , and 
introduce 

Z = -^JDJ = (JS)(SJ) T (5) 

where Ski is the ^-th component of Sk and J = (Jki) = 
{5kr - 1 IN) is called the centering matrix whose effect for S 
is (JS)ke — Ski ~ YikSktlN. Next we decompose Z into its 
eigenvalues and eigenvectors as Zp l = Ajp t with A t > A i+l . 
Then we get Z = PAP T where A;; = AjSij and Pn is the 7-th 
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FIG. 5. Neighboring graph Q. (a) Schematic illustration of the 
sample points (o) and the 3-neighboring graph Q (solid line). The 
heavy solid line is the shortest path between points i and j and the 
dashed line is the Euclidean distance, (b) 500 sample points (+) and 
4-neighboring graph Q (line) for v = 0.1130. (c) Graph Q for 2- 
band chaos (vp =0.1128). The two clusters Bi, B2 caused by 2-band 
chaos are connected by the shortest path C. 



component of p r Therefore, we obtain the matrix 
U = PA 112 , (A 1/2 ), v = V^0' 



(6) 



which corresponds to the matrix JS and the new point «, 
which is the z'-th row vector of U. Clearly the new points 
h, are reconstructions of the original points s, and recover the 
distance \Sj - s } \. If A„ » /t„+i, we can approximate U by 
its projection into the subspace spanned by the eigenvectors 
{p l ,p 2 ,...,p n }. 

As we can find the points «, only from the distances, we 
start with the geodesic distances dtj instead of |s; - Sj\ and 
get the new low dimensional vector m, corresponding to the 
z'-th sample point on the filament. Finally, we obtain n- 
dimensional representations «, = (u^\ u. , . . . , u.) whose 
distances preserve the geodesic distances in the original space 
as much as possible. This procedure is essentially the same as 
the principal component analysis for the data matrix U 0>0]. 



B. Results for the mass-spring model 

Figure [6] (a) shows configurations of sample points for the 
dripping faucet data (t„_i,t„) onto the (m (1) ,m (2) ) plane ob- 
tained from ISOMAP. Here, the spread of points in the direc- 
tion of the m (2) component is much smaller than that of the u m 
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FIG. 6. Results of the application of ISOMAP for Fig. [5](c). (a) 
Configurations of sample points onto the (u {l \ i( <2) ) plane, (b) Plot of 

sample points in the (t b _i,t„, m^'J space and its projection onto the 
(t„_i,t„) plane, (c) Relation between the new coordinate u and the 
droplet mass m. 



component (about 0.6%). This means that the points on the fil- 
ament are almost explained only by the first component m (1) , 
which implies that ISOMAP succeeds to unfold the attractor 
to a straight line and u m is considered as a new coordinate 
along the filament (Fig. |6|b)). Hereafter, u„ is abbreviated 
to u„ as we mainly use the first component. The successful 
unfolding can be also confirmed from the one-to-one relation 
between the mass m„ and the new coordinate u„ in Fig. |6fc). 

As any point on the filament is deterministically mapped to 
another point, the time evolution of u is described by a one- 
dimensional map u„ = g{u„-\) which is shown in Fig|7j In 
addition, the points whose spread in the m (2) direction in Fig. 
|6ja) is relatively large are located around the folding point 
(critical point) of the one-dimensional map in Fig. [3a). 

The results in Fig|7]show the expected relationship between 
m„_i and m„. We approximate this one-dimensional map g(u) 
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FIG. 7. (a),(b) Return plots of the new variable m (i) generated by 
ISOMAP for Fig.|5](b) and (c), respectively. 



by a locally quadratic function 



(7) 



where a/u) is determined by the locally least square method, 
i.e. aj(u) satisfies 



min Y/g(M«-i)-M«) 2 expj-(— — -) 



and we use cr = 0.05. 



C. Statistical properties 

From the one-dimensional map #(m), we can get the natural 
invariant density p(u) which is the base for the discussion of 
the statistical properties. Here, p(u) satisfies the equation 



p(u) = H(p(uJ) = J p(v)6(u - g(v))dv 

- Z 



u„:u=g(u„) 



p(u n ) Ag 

, g (u) — — 

\g'(u n y\ du 



(8) 



where H is called Frobenius-Perron operator and p(u) 
is generally expected to be the empirical distribution 
limA^oo 2„=i <5( M ~~ u„)/N for a chaotic orbit. 

To solve approximately Eq.®, we divide the domain of 
g(u) into the intervals I,, (/ = 1,2, . . .) whose edges are in- 
verse mapping points g~ n (u*),(n = 0, 1,2, . . .), where u* is 
the critical point (the minimum in Fig|7](a) and the maximum 
in Fig|7](b) ) of the function g(u) and expand p(u) as 



PO) = J] aidiu), ei(u) = j* 
We substitute Eq.® in Eq.© and get 
a; = YH u aj, Hjj = 'j- 

V \g'U 



1,(m € L) 
(others). 



(9) 



1, if I y c gild 



| '(»f)l' Aj \0,iflj£g(ld (10) 

where u\ is the center of Iy. The solution is given by the 
eigenvector corresponding to the eigenvalue 1 of the matrix 
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FIG. 8. (a) Relation between the new valuable u and the dripping 
interval r. (b) Distribution function of dripping intervals P(t) which 
is derived from Eq. lll2t . The result from the direct simulation (10 6 
droplets) are also plotted by the symbol +. 



H. The above-mentioned method is a kind of the Galerkin- 
approximation which is often used in the study of the one- 
dimensional map IllOll . 

The internal state variable u is not a natural physical quan- 
tity for the dripping faucet system. However we can derive 
any physical quantity A from m, because the quantity A on the 
filament is determined by m, i.e., A = A(u) and its long-time 
average is calculated by (A) = J A(u)p(u)du. Actually, as the 
results of ISOMAP provide the relationship between u n and 
(t„_i, t„) as shown in FigJHIa), the important observable vari- 
able t of the dripping faucet system is determined by 



t = (f>(u) 



(11) 



where the function <f)(u) is approximated in the same way as 
g(u) (see Eq.(|7]i). First we get the distribution function 



P(T) = 



P(u) 

I0'(")l 



(12) 



which is one of the most basic properties of the dripping 
faucet. The result in Fig. [8]shows P(t) with many peaks which 
are generated by the folding processes of g(u). These proper- 
ties are consistent with the result from the direct simulation of 
Eq.© and ©. 

Next, we calculate the Lyapunov exponent Hill and 
the topological entropy lfl2ll which characterize the 
stability and the variety or complexity of chaotic or- 
bits, respectively. The Lyapunov exponent of g(u) 
is defined by A = limAr^o^l/AOlnldiWdMol = 
liniiv^ooO/AO Z*=0 ln|dM„+i/d«„| and given by 



(ln|^'(«)|). 



(13) 



The topological entropy is equal to the largest eigenvalue of 
the transfer matrix (3 in Eq.(fT0l>. As both quantities are invari- 
ant under the transformation from u to m, the results from g(u) 
should coincide with the Lyapunov exponent and topological 
entropy of f{m) (if the map g{u) is an appropriate descrip- 
tion of the dripping faucet dynamics given by Eqs. (0 and 
©). They are cited in Table J] The Lyapunov exponents are 
in good agreement, but the topological entropy from g(u) is 
slightly smaller than the results from f(m). This may mean 
that the sample points t„ do not include rare orbits, because 
the results are based on only 500 sample points. 
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TABLE I. Lyapunov exponents and topological entropies calculated 
from f(m) and g(u). 



vo 


f(m) or g{u) Lyapunov Exp. topologcal entropy 


0.1128 


/(«) 


0.253 


0.859 




g(») 


0.253 


0.855 


0.1129 


f{m) 


0.306 


0.913 




g(») 


0.291 


0.882 


0.1130 


f(m) 


0.350 


0.948 




g(») 


0.346 


0.903 



0.5 



C(n)/C*(0) o 



-0.5 



10 20 30 40 50 60 



n 



shown in Fig. [TOfc). The absolute value of its slope is nearly 
one in almost all regions i.e. the instability of the orbits is 
weak which is related to the fact that the time series contains 
one or two periodic like motions. 

It has been pointed out that ISOMAP is topologically un- 
stable for small noise [13]. In actually, the neighboring graph 
Q around the folding point of the filament is affected by ex- 
perimental noise or high-dimensional dynamics and has some 
short cuts which are out of the filament. Therefore, the first 
return map of u has a muliti-valued structure around the crit- 
ical point (maximum point) of g{u). However, this effect is 
small and localized so we can say our method is a promising 
method for not only the numerical study but also experimental 
data. 

Our experiment is a first step and currently we cannot keep 
it stationary to measure long time series. We work on a revised 
set-up and we shall present extended experimental results in a 
future article. 



IV. SUMMARY AND DISCUSSION 



FIG. 9. Time correlation function of r„. The results from g{u) (o) 
and the direct simulation (10 6 droplets) from Eq. ((3} and ©(x) are 
plotted. The two curves show the decay for a period-2 and a long 
period (about 23) oscillation. 



We can calculate the time series r„ from the one- 
dimensional map (0 and Eq.dTTb. but its long time behavior 
has a large difference from the direct simulation of Eq.(|3) and 
©, because the dynamics is chaos. We also calculate the time 
correlation function of t„ 




= <0(« W(«))> - (<K") 2 ) (14) 

where f = Urrijv-^l/jV) £m=o T m . The result in Fig|9]shows 
the exponential decay for a period-2 and long period (about 
23) oscillation. Their properties are also shown in the result 
which is calculated directly from t„ by the simulation of Eq. 
(0 and |@). The good coincidence shows that by using the in- 
ternal variable u, we can discuss not only the static properties 
but also the dynamical property of the original observable r. 

D. Application to real experimental data 

Last we show preliminary results of the application of the 
dimension reduction method to our real experimental data. A 
short time series whose water flux is almost stationary and its 
first-return plot r„_i vs. t„ are shown in Figs. [Tol a) and (b), 
respectively. This return plot shows a multi-valued function 
in the wide region and cannot describe the dripping faucet dy- 
namics. The application of our method leads to the internal 
variable u and the one-dimensional map u„ = gf(«„-i) which is 



In this paper, employing the dripping faucet system as an 
illustrative example, we studied the problem of constructing a 
surrogate variable for the internal state of the chaotic dynami- 
cal systems from time series using manifold learning analysis. 
Especially, when the time-delay embedding of the observed 
time series forms a one-dimensional curved structure, we suc- 
ceeded to obtain one-dimensional deterministic maps associ- 
ated with the surrogate variable. The statistical properties of 
the original chaotic system were successfully reproduced by 
its surrogate system. 

In real-world applications, not all original state variables of 
the system can be directly observed. Instead, only some of the 
original state variables or their transformations are observed in 
experiments. So, it is often seen that the manifolds obtained 
from the time-delay embedding may be very complex even if 
their dimensionality is low, which leads to multi-valuedness in 
the return plots [5]. For example, the time series of inter-spike 
intervals is mainly observed in neural systems and its return 
plot often exhibits a multi-valued function lfl4ll . Besides neu- 
ral systems, there is a number of such examples, e.g., laser 
systems 11151 llq] , passive biped walkers [17], and social ac- 
tivity models 111 811 . Extracting the deterministic relationship 
from the observations of these models may be also done using 
dimension reduction methods. 

Besides, the return plot of u n -\ vs. u„ (Fig. |7) is less 
smooth compare to that of m n _i vs. m„ (Fig. [3] (a)). This is 
because methods of manifold learning are generally unsuper- 
vised ones, just using the information of t„. If the assump- 
tion that the training data is generated from a dynamical sys- 
tem with a simple mapping form say the logistic parabola can 
be incorporated additionally to the manifold learning as some 
constraint or penalty terms, we can obtain a more refined re- 
turn plot which may be more interpretable to us. We would 
like to develop such problems in our future works. 

In this paper, we have been only concerned with the case in 
which the internal state behind the observed time series obeys 
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(a) 



0.32 
0.31 
0.3 
0.29 



0.32 



0.31 



+ Vv^^ViM^^-t. + + + + ++++ + 4 



0.29 



50 100 150 200 

n 




0.3 0.31 
Tn-1 



0.32 



a one-dimensional dynamical system, this is of course an ideal 
case. We have to extend our approach to higher dimensional 
cases (two, three dimensional maps). As formal method- 
ologies (application of ISOMAP or other manifold learning 
methods) are not restricted to one dimensional manifolds, the 
presented approach can in principle be extended to higher di- 
mensional cases. It should be noted, however, that ways of 
acquiring training samples to obtain a lower dimensional rep- 
resentation become more important. For example, let us con- 
sider the situation where the Henon attractor J[ is lying on a 
two dimensional nonlinear manifold M embedded in, say, the 
three dimensional Euclidean space R 3 . In order to obtain a 
lower dimensional representation of At, not only the data on 
3K, but also the data associated with transient dynamics are 
needed because £R is too thin to recover the whole two dimen- 
sional structure of M. In addition, the non-uniformity in the 
natural measure on J?l and its transient area affects the perfor- 
mance of manifold learning. 



0.04 
0.02 
U n 
-0.02 
-0.04 



X 



-0.04 -0.02 0.02 0.04 
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FIG. 10. Results for our experimental data, (a) Dripping time series 
t„ which contains one or two periodic motion around n = 50, 170. 
(b) First-return plot of r. It shows a multi- valued function in the right 
region, (c) First-return plot of inner state variable u which suggests 
that there is a one-dimensional map g{u). It may has a wavelike 
pattern in the right region. 



ACKNOWLEDGMENTS 



This study is partially supported by Grant-in-Aid for Sci- 
entific Research (No. 22740258), the Ministry of Education, 
Science, Sports, and Culture of Japan. The research leading 
to the results has received funding from the European Com- 
munity's Seventh Framework Programme FP7/2007-20 13 un- 
der grant agreement No. HEALTH-F2-2009-241526, EU- 
TrigTreat. Furthermore, support by the Bernstein Center 
for Computational Neuroscience II Gottingen (BCCN grant 
01GQ1005A, project Dl) is acknowleged. H.S. is grateful to 
S. Akaho for fruitful discussions and comments. 



[1] R. Shaw, The Dripping Faucet as a Model of Chaotic System 

(Aerial Press, Santa Cruz, 1984). 
[2] K. Kiyono and N. Fuchikami, J. Phys. Soc. Jpn. 68, 3259 

(1999). 

[3] N. Fuchikami, S. Ishioka, and K. Kiyono, J. Phys. Soc. Jpn. 68, 
1185 (1999). 



[4] T. Katsuyama and K. Nagata, J. Phys. Soc. Jpn. 68, 396 (1999). 
[5] F. Christiansen, P. Cvitanovic and V. Putkaradze, Nonlinearity 
10, 55 (1997). 

[6] T. Hastie, R. Tibshirani, J. Friedman, The Elements of Sta- 
tistical Learning: Data Mining, Inference, and Prediction 
(Springer- Velag, New York, 2001). 



9 



[7] J.B. Tenenbaum, V. de Silva and J.C. Langford, Science 290, 
2319 (2000). 

[8] T.F. Cox and M. A.A. Cox, Multidimensional Scaling (Chapman 

and Hall, London, 2000). 
[9] C. Williams, Machine Learning 46 11 (2002). 
[10] T. Kohda and K. Murao, IEICE E73 793(1990). Erik Bollt, 

Pawel Gora, Andrzej Ostruszka and Karol Zyczkowski, SIAM 

J. on Applied Dynamical Systems, 7, 341(2008). 
[11] E. Ott, Chaos in Dynamical Systems (Cambridge University 

Press, Cambridge, 2002). 
[12] B.-L. Hao, W.-M. Zheng Applied Symbolic Dynamics and 

Chaos(World Scientific Publishing Company, Singapore, 

1998). 



[13] M. Balasubramanian and E. L. Schwartz, Science 295, 7, 
(2002). 

[14] U. Feudel, A. Neiman, X. Pei, W. Wojtenek, H. Braun, M. Hu- 
ber and F. Moss, Chaos 10, 231 (2000). 

[15] U. Hubner, C.-O. Weiss, N.B. Abraham and D. Tang, In: 
A.S. Weigend and N.A. Gershenfeld (eds.), Time-Series Pre- 
diction: Forecasting the Future and Understanding the Past, 
73, (Westview Press, Boulder, 1993). 

[16] J. Used and J.C. Martin, Phys. Rev. E 79, 046213 (2009); ibid. 
82, 016218 (2010). 

[17] A. Goswami, B. Thuilot and B. Espiau, Int. J. of Robotics Res. 
17, 1282 (1998). 

[18] G. Feichtinger, L.L. Ghezzi and C. Piccardi, Int. J. Bifurcation 
and Chaos 5, 255 (1995). 



